#==========================================================
# Simulation of an AR(1) Model using arima.sim()
# Code is from Z&W (slightly modified)
# Also illustrates par(mfrow(2,2)) and graph labeling 
#=======================================================


#The simulation using start.innov burn-in. We simply add "something" (here it is 1) 
#to get a series that has a mean differing from zero.


e<-rnorm(100,sd=1)
e.start=rnorm(25, sd=1)
y.ar1<- 1+arima.sim(model=list(ar=0.75), n=100,innov=e, start.innov=e.start)

#Here the values in start.innov are used for the "burn in." For a nice AR(1) model
#a burn-in period of 25 probably suffices, but for more complicated models you may want a 
#longer burn-in period.

#By the way, you "need" just as much burn-in to generate a series of length 100 as to 
#generate one of length 5000. 

# Let's check on the mean to see if it makes sense. 
# We may run our simulation several times.

mean(y.ar1)

# If we want to have the simulation begin with the same random seed each time
# we can explicitly set the pseudo-random number generator seed

set.seed(101)

# When we use this seed we get a larger mean than we might have expected, but this
# but due to the reasonably large variance of this series (calculate tau!), the mean
# will bounce a lot from seed to seed. 

#Now lets do some pictures --- and put three on just one page:
#
#(1) the simulated series,
#(2) the theoretical ACF
#(3) the empirical ACF for our simulation

par(mfrow=c(2,2))
tsplot(y.ar1, main="Simulated AR(1)")
abline(h=1)

gamma.j<-rep(.75,10)^(0:9)
tsplot(gamma.j, type="h", main="Theoretical ACF for AR(1)", ylab="Autocorreleation", xlab="Lag+1")
temp<-acf(y.ar1, lag.max=10, plot=F)
tsplot(temp$acf,type="h", main="Empirical ACF for Simulated AR(1)", 
ylab="Autocorreleation", xlab="Lag+1")
 

#One line to make it prettier:
abline(h=0)

#How do the series compare? To me the typical run looks pretty reasonable, 
#though you may notice that at higher lags the sample ACF tends to "turn up."
#Thus, there is a sampling effect that is not present in the theoretical ACF

#You might want to look at the help screens for the acf function.